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20.  cont'd 

radiation  model  has  been  incorporated  into  the  one-dimensional 
model  to  permit  the  radiative  flux  divergence  term  to  be 
computed  in  a fully  coupled  manner.  Sample  calculations  have 
been  made  both  to  verify  that  the  program  behaves  in  a reason- 
ably physical  manner  and  to  exemplify  some  of  the  types  of 
boundary  layer  distributions  which  may  be  expected  to  occur. 

The  coupled  thermal  radiation  model  was  used  to  make  calcula- 
tions of  the  typical  diurnal  variations  in  the  boundary  layer 
over  both  land  and  water.  The  time  variation  of  the  boundary 
layer's  stability  is  quite  different  over  water  than  it  is  over 
land.  This  leads  to  sea-breeze  circulations  in  the  coastal 
planetary  boundary  layer  which  are  exemplified  in  calculations 
using  the  two-dimensional,  unsteady  version  of  the  model. 
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I . INTRODUCTION 


This  report  describes  the  progress  made  over  the  past  year 
toward  the  development  of  a program  for  computing  the  detailed 
low-level  atmospheric  distributions  of  velocity,  temperature, 
moisture,  refractive  index,  and  the  turbulent  variances  of  these 
quantities  for  marine  environments.  The  program  is  designed  to 
provide  the  capability  of  taking  the  predicted  output  of  a large- 
scale  meteorological  forecast  model,  such  as  that  developed  at 
Fleet  Numerical  Weather  Central,  and  using  this  as  boundary  condi- 
tions to  make  predictions  of  the  detailed  microstructure  of  the 
planetary  boundary  layer  at  desired  locations.  This  program 
development  also  provides  a means  for  checking  and  upgrading  the 
boundary  layer  parameterization  used  in  the  global  forecast  model 
which  should  ultimately  permit  ai.  improved  accuracy  in  the  large- 
scale  forecast. 

Our  approach  to  this  program  has  been  to  use  the  invariant 
second-order  closure  model  of  turbulence  developed  by  Dr.  Coleman 
duP.  Donaldson  and  his  associates  at  A.R.A.P.  over  the  past  few 
years.  The  fundamentals  of  this  approach  are  given  in  Ref.  1,  A 
review  of  the  status  of  this  model  as  applied  to  a v;ide  variety  of 
turbulent  flow  problems  is  given  in  Part  II  of  Ref.  2.  Particular 
applications  of  the  model  as  applied  to  atmospheric  problems, 
including  comparisons  with  experimental  data,  are  documented  in 
Refs.  2-6. 

Reference  7 is  a technical  report  detailing  the  model  develop- 
ment, sample  calculations,  and  verification  comparisons  made  under 
our  initial  contract.  It  describes  the  addition  of  humidity  and 
the  second-order  turbulent  correlations  Involving  humidity  as 
variables  to  our  dry  atmospheric  boundary  layer  program  previously 
deve^'ped  for  the  Environmental  Protection  Agency.  These  variables 
were  ‘ ided  v/ithout  the  need  to  add  any  new  coefficients  to  the 
model  by  observing  that  in  the  limit  of  vanishingly  small  fluctua- 
tions in  either  humidity  or  temperature  both  are  transported  as 
passive  quantities  and  thus  should  be  modeled  similarly.  The 
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different  influence  of  humidity  and  temperature  on  buoyant  pro- 
duction of  turbulence  is  Included  explicitly  in  the  equations 
fur  the  turbulence,  with  no  need  to  model  these  terms.  Verifi- 
cation o»'  the  model  predictions  for  humidity  transport  were  made 
by  compai'ison  with  evaporation  experiments  in  a wind  tunnel 
(Ref,  8). 

Tht  boundai-y  conditions  at  the  surface  of  the  earth  were 
modified  to  provide  for  appropriate  coupling  to  the  sea  state 
by  using  Froude  number  scaling  following  Wu  (Ref  9)  to  provide 
a relationship  between  the  effective  aerodynamic  roughness  of 
the  sea  surface  and  the  surface  shear  stress  velocity,  This 
allows  both  the  low-level  wind  distribution  and  the  surface 
roughness  to  be  determined  by  the  model. 

The  total  water  content,  liquid  plus  gas,  is  transported  by 
the  turbulence.  Provision  for  the  water  content  to  change  its 
phase  has  been  made  by  assuming  that  when  both  liquid  and  water 
vapor  coexist  their  relative  proportion  is  determined  by  equili- 
brium saturation  conditions.  The  heat  released  in  the  conden*r 
sation  process  is  included  an  both  the  mean  average  temperature 
equation  and  in  the  correlation  equations  Involving  temperature 
fluctuations . 

Using  the  predicted  distributions  of  temperature,  humidity, 
and  pressure,  a calculation  of  the  modified  refractive  index,  M , 
has  been  incorporated  in  the  program.  Local  minlmums  in  the  M 
distribution  with  respect  to  altitude  directly  indicate  the 
presence  of  a radar  duct.  Since  we  are  predicting  the  second- 
order  correlations  oetv/een  the  turbulent  fluctuations  in  tempera- 
ture and  humidity  as  well  as  the  average  scale  of  the  turbulence, 
we  have  available  the  information  to  also  compute  the  structure 
of  the  fluctuations  in  refractivity , 

Reference  7 includes  the  results  of  several  sample  calcu- 
lations; e.g.:  (a)  a sample  calculation  using  output  from  PNWC 
(supplied  by  J,  Kaitala)  as  upper  and  surface  boundary  conditions 
for  our  boundary  layer  program;  (b)  a comparison  of  predicted 
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temperature  structure  parameters  with  the  observations  of 
Wyngaard,  et  al,  (Ref.  10);  (c)  a calculation  with  boundary 
conditions  roughly  corresponding  to  the  conditions  observed  in 
the  Atlantic  Tradewind  Experiment  by  Augstein,  Schmidt,  and 
Ostapoff  (Ref  11);  and  (d)  a calculation  simulating  shoreline 
conditions  ior  either  a dry  land  breeze  over  the  sea  or  a moist 
sea  breeze  over  dry  land. 

This  report  describes  two  major  modifications  to  the  model 
described  in  Ref.  7.  These  are:  (a)  the  1 .crease  in  the  dimensions 
of  the  program  to  a two-dimensional,  unsteady  calculation  to  per- 
mit the  prediction  of  shoreline  conditions  developing  in  time,  and 
(b)  the  incorporation  of  the  radiative  flux  divergence  term  into 
the  one-dimensional  system  of  equations  in  a coupled  manner.  The 
numerical  developments  required  to  extend  the  program  to  two  di- 
mensions are  discussed  in  the  next  section. 

The  radiation  model  coupled  into  the  one-dimensional,  boundar*' 
layer  model  is  described  in  Appendix  A.  The  incorporation  of  a 
coupled  radiation  model  is  important  under  stable  atmospheric  con- 
ditions when  the  comparative  ratio  of  the  divergence  of  radiation 
heat  flux  to  that  of  the  turbulent  heat  flux  may  reach  order  one. 
The  primary  coupling  between  the  turbulent  transport  and  radiation 
comes  through  the  humidity  distribution.  The  water  vapor  content 
has  a strong  influence  on  the  long-wave  radiative  cooling,  while 
the  liquid  water  content  is  the  most  dominating  factor  in  the 
short-wave  radiative  heating.  Verification  comparisons  with 
observational  data  are  included  in  Appendix  A, 

We  have  also  made  some  modifications  to  the  condensation 
algorithm  to  provide  for  a transition  regime  betv^een  the  completely 
unsaturated  and  completely  saturated  regimes.  The  need  for  and 
consequences  of  this  modification  are  discussed.  Calculations  have 
been  made  both  to  verify  that  the  modified  program  behaves  in  a 
reasonable  fashion  and  to  exemplify  some  of  the  types  of  boundary 
layer  distributions  which  may  be  expected  to  occur. 
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Sample  caiculaticns  of  -e  diurnal  variation  in  the  planetary 
. -undary  layer  over  Doth  land  and  water  have  been  calculated  wic> 
detailed  thermal  radiation  coupling.  The  conditions  for  the 
calculation  over  land  are  picked  to  correspond  to  that  for  a 
typical  midv;estern  summer  day  previously  calculated  ignoring 
radiation  coupling  and  published  in  Refs.  2 and  This  permits 
the  significance  of  the  radiation  coupling  for  such  a calculation 
to  be  directly  assessed,  A sample  calculation  over  a boundary 
layer  over  water  is  then  presented  to  show  the  interaction  be- 
tween radiation,  cloud  dynamics,  and  turbulence.  These  two  sample 
diurnal  variation  calculations  clearly  show  that  the  time  variation 
of  the  boundary  layer’s  stability  Is  quite  different  over  water 
than  it  is  over  land. 

Finally,  the  two-dimensional,  unsteady  version  of  the  model 
has  been  used  to  calculate  the  typical  variation  in  the  coastal 
planetary  boundary  layer.  The  results  of  this  calculation  have 
been  written  for  presentation  at  the  AMS  Conference  on  Coastal 
Meteorology  to  be  held  at  Virginia  Beach,  September  21-23>  1976. 
This  paper  is  included  as  Appendix  B to  this  report.  The  re- 
sulting diurnal  variation  in  the  sea  breeze  induced  by  the  strong 
stability  differences  In  the  boundary  layer  response  over  the  land 
and  that  over  the  water  produces  a strong  asymmetry  between  the 
sea-breeze  and  the  land-breeze  circulation  patterns.  In  previous 
sea-breeze  models  it  was  necessary  to  impose  eddy  viscosities 
which  were  a strong  function  of  time  and  space  to  gain  this 
asymmetry,  In  the  present  model  it  is  achieved  without  the  need 
tc  introduce  any  new  empirical  information. 


II.  MODEL  EXTENSION 


Two  major  modifications  to  the  model  have  been  made.  The 
numerics  of  the  model  have  been  extended  to  permit  two-dimensional, 
unsteady  calculations,  and  a radiative  flux  divergence  calculation 
has  been  added  to  the  one-dimensional  system  of  equations  to  permit 
a direct  coupling  between  the  distributions  of  temperature  and 
water  vapor  predicted  by  turbulent  transport  and  the  radiative 
energy  source. 

The  only  other  change  to  the  model  during  the  past  year  has 
been  a refinement  in  our  condensation  algorithm  to  provide  for  a 
transition  regime  between  the  completely  unsaturated  and  completely 
saturated  regimes.  Sample  calculations  had  shown  that  this  was 
desirable  to  simulate  conditions  corresponding  to  scattered  cloud- 
iness, The  new  algorithm  permits  condensation  to  be  Influenced  by 

the  turbulent  fluctuations  in  H and  H„  as  well  as  by  their 

s 

ensemble-average  value. 

The  development  and  verification  of  the  radiation  model  are 
given  in  Appendix  A, 

2.1  Program  Extension  to  Two  Dimensions 

To  include  the  physics  generated  by  the  true  ellliptic  nature 
of  the  shoreline  sea-breeze  problem,  we  must  allow  additional 
variation  in  the  direction  normal  to  the  shore.  The  program  used 
to  develop  the  solutions  discussed  later  and  in  Appendix  B is  an 
extension  of  a computer  program  developed  with  the  support  of  EPA 
through  our  continuing  turbulent  modeling  effort.  This  program 
initially  solved  for  the  primitive  variables  (W,  H,  6^,  U,  V)  and 
the  turbulence  in  the  y-z  computational  mesh  in  a time-dependent 
mode.  In  this  scheme,  the  vertical  velocity  W was  obtained  using 
the  boundary  layer  approximations  so  that  it  was  solved  by  quadra- 
ture from  the  horizontal  velocity  V through  continuity.  A review 
of  the  solution  process  may  be  found  in  Ref.  2. 

To  complete  the  sea-breeze  circulation,  it  is  necessary  to 
provide  a consistent  V/  boundary  condition  at  the  top  of  the 


slmalation  domain.  The  aoe  of  the  continuity  equation  to  deter- 
mine W allows  the  use  of  only  one  boundary  condition  (namely, 
v;  = 0 at  z = 0)  and  forces  VJ  to  be  nonzero  at  the  top  of  the 
computational  domain.  To  permit  the  imposition  of  W = 0 at  the 
top  (and  thereby  closing  the  domain),  we  developed  a hybrid  solu- 
tion scheme  involving  the  streamfunction  . The  horizontal 
velocity  V is  solved  at  each  time  step  by  the  horizontal  momentum 
equation.  The  vorticity  is  found  by  differentiating  V in  the 
vertical  direction;  the  streamfunction  i{)  is  obtained  by  a direct 
solution  of  Poisson's  equation  using  the  solver  developed  by  Swarz- 
trauber  and  Sweet  (Ref.  12).  The  top  boundary  value  of  ij/  is  taken 
as  constant,  forcing  V/  to  be  zero  along  the  top  boundary.  The 
vertical  momentum  equation  provides  the  equation  for  the  hydrostatic 
pressure.  The  resulting  horizontal  pressure  gradients  are  the 
forcing  functions  for  the  sea-breeze  circulation. 

The  greatest  difficulty  in  the  numerics  of  this  two-dimensional, 

unsteady  system  is  associated  with  the  internal  waves  in  the  region 

of  the  stable  temperature  inversion.  Since  these  stable  potential 

temperature  gradients  can  become  relatively  strong  at  the  top  of  the 

mixed  layer,  the  characteristic  Brunt-Vaisala  period  can  be  of  the 

order  of  tens  of  seconds.  Since  the  full  turbulent  correlation 

equations  follov/  these  waves,  it  is  necessary  to  run  with  very  small 

time  increments.  This  makes  a full  day's  simulation  of  the  two- 

dimensional  coastal  boundary  layer  unduly  lengthy.  To  circumvent 

this  problem,  the  sample  calculations  shown  here  and  in  Appendix  B 

have  been  run  with  upwind-differencing  and  with  quasi-equilibrium 

turbulence.  Upwind-differencing  Introduces  some  numerical  damping 

into  the  calculation  so  that  larger  time  steps  may  be  taken.  The 

quasi-equilibrium  turbulence  approximation  is  detailed  in  Ref.  2. 

The  basic  approximation  is  to  drop  out  all  time  derivative  and 

diffusion  terms  in  the  second-order  correlation  equations,  except 
2 

in  the  q equation.  This  has  the  desirable  effect  of  filtering 
out  the  wave  dynamics  from  the  correlation  equations.  The  penalty 
for  this  is  some  loss  in  the  faithfulness  with  which  the  equations 
model  the  turbulence. 
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We  cannot  give  a precise  bound  on  the  errors  Introduced  into 
the  calculation  by  these  approximations  but  believe  that,  at  least, 
the  qualitative  behavior  of  the  results  is  correct.  Definitive 
accuracy  estimates  must  await  comparison  with  reliable  field  data. 

It  should  perhaps  be  restated  here  that  the  two-dimensional 
program  does  not  currently  have  either  the  coupled  radiation 
aalculitlon  or  the  condensation  algorithm  Included.  These  two 
physical  phenomena  are  being  checked  out  in  the  simpler  one- 
dimensional  program  before  inclusion  into  the  two-dimensional 
system. 

Figure  1 demonstrat.  the  difference  between  the  results  of 
a two-dimensional,  elliptic  calculation  for  the  shoreline  boundary 
layer  and  that  of  a two-dimp  'sional,  parabolic  calculation  in  even 
the  steady  case.  For  this  calculation,  the  geostrophic  wind  is 
assumed  to  be  directed  normal  to  the  shore  with  the  water  temper- 
ature 11®C  cooler  than  the  land.  The  inflow  conditions  are  held 
constant  at  conditions  corresponding  to  a typical  quasi-steady 
boundary  layer  over  water  with  an  inversion  height  of  approximately 

1500  m. 

This  problem  may  be  solved  using  the  one-dime  onal  program 
by  assuming  that  the  flow  is  parabolic  and  marchinti  in  the 
direction  normal  to  the  shoreline.  The  contours  of  constant  turbu- 
lent kinetic  energy  obtained  in  this  way  are  shown  in  Fig.  1-a. 

This  calculation  does  not  permit  any  Influence  of  recirculations. 
Figure  1-b  shows  the  corresponding  variable  as  predicted  by  the 
two-dimensional  elliptic  calculation.  This  is  obtained  by  starting 
the  two-dimensional  calculation  with  the  conditions  obtained  from 
the  parabolic  calculation  and  then  marching  in  time  until  a type  of 
steady  state  is  reached.  Figure  1-b  shows  more  structure  than  can 
be  observed  In  Pig.  1-a.  Figure  1-c  is  the  qualitative  field 
observations  as  obtained  from  multiple  passes  with  an  aircraft 
(Ref.  13)  along  the  shore  of  Lake  Michigan.  The  value  of  the 
geostrophic  wind  at  the  time  of  their  flight  is  not  specified,  so 
it  is  not  possible  to  make  a quantitative  comparison.  The  height 


Figure  1-b . Boundary  layer  growth  in  turbulent  kinetic  energy  along  the  shoreline  with 
the  geostrophic  wind  from  the  water  at  10  m/sec  and  the  land  warmer  than 
the  water;  elliptic  calculation. 


01  the  boundary  layer  thickness  at  50  km  inland  would  increase 
the  geuotx’ophic  wind  uin  ^ decreased  for  a specified  temperature 
difference  between  the  land  and  the  water.  We  did  use  the  same 
temperature  difference  in  the  calculation  as  was  specified  in  the 
ooservations. 

Both  the  elliptic  predictions  and  the  observations  show,  at 
least  qualitatively,  the  same  structure  to  the  turbulent  kinetic 
energy  field.  The  local  maximum  in  nearest  to  the  shoreline 

is  predicted  by  both  the  elliptic  and  the  parabolic  calculation. 

This  increase  in  turbulence  is  caused  by  the  increased  shear  stress 
and  heat  flux  in  the  internal  boundary  layer  which  begins  at  the 
shoreline.  As  the  new  boundary  layer  increases  in  thickness,  the 
shear  decreases  and,  consequently,  so  does  the  turbulent  kinetic 
energy.  The  added  detail  of  the  local  recirculations  appears  to  be 

p 

necessary  to  predict  tne  second  local  maximum  in  q , The  vertical 
potential  temperature  contours  for  this  flow  are  shovm  in  Pig.  2. 
Again,  the  undulating  character  of  the  boundary  layer  growth  is 
clearly  in  evidence.  As  previously  explained,  the  results  of  this 
computation  cannot  be  used  to  quantitatively  verify  the  program, 
but  it  does  show  that  the  results  are  at  least  qualitatively  correct. 

2.2  Turbulent  Contribution  to  Condensation 

As  detailed  in  Ref.  7,  the  change  of  phase  of  water  was 
originally  taken  as  entirely  governed  by  the  mean  values  of  H and 
H . As  explained  there,  we  were  relying  on  turbulent  diffusion 

O 

in  the  equations  to  smooth  out  the  transition  between  saturated 
and  unsaturated  conditions.  Under  conditions  when  light, scattered 
clouds  might  be  anticipated,  this  appeared  to  cause  some  problems. 

The  program  attempted  to  compensate  for  its  inability  to  specify 
spatiaily  scattered  clouds  by  prescribing  a birth  and  decay  of 
clouds  in  time  over  a period  of  the  order  of  minutes. 

Under  the  assumption  of  equilibrium  condensation,  the  change 
of  phase  should  be  governed  by  the  instantaneous  values  of  H and 
H_  . Since  liquid  water  content  ? is  physically  a positive 

O 
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definite  quantity,  it  may  be  seen  that  the  contributions  made 

when  H + H'  > H + H’  cannot  be  balanced  by  contributions  from 

H + H’  < H + H'  . Thus, 
s s ’ 

C = H - Hg  (1) 

only  when  the  flow  Is  sufficiently  saturated  that  no  local  unsatu- 
x-ated  region  is  contained  within  an  eddy. 

To  approximately  account  for  this  phenomenon,  we  have  added  a 
transition  regime  between  the  completely  unsaturated  and  completely 
saturated  regimes.  This  occurs  when  the  fluctuations  in  H and 
may  be  large  enough  to  make  the  sign  of  H - to  the  sign 

S _ ^ o 

of  H - H . letting  the  fluctuations  be  characterized  by  their 
s 

rms  values,  the  transition  regime  may  be  defined  as  occurring 
whenever 


To  precisely  specify  the  liquid  water  content  in  this  regime,  it 
would  be  necessary  to  k.iow  the  distribution  of  H'  and  H‘  about 

. , o 

H and  . Since  this  is  more  information  than  is  possible  in 
a model  based  on  ensemble  averages,  we  will  approximate  the  average 
liquid  water  content  in  this  region  as 


5 


(3) 


Equation  (3)  has  the  correct  limits  of  zero  and  H - H when  the 
equal  sign  in  Eq.  (2)  holds.  Further,  when  H = H and  H*  (or 

D 

H’)  is  zero,  it  recognizes  that  only  negative  H'  (or  positive 
s s 

H')  will  lead  to  condensation. 

The  coefficients  , related  to  the  effective  increase  in 
the  specific  heat  capacity  of  the  mixture  due  to  change  of  phase, 
and  K2  , the  effective  change  in  the  neutral  lapse  rate  due  to 
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change  of  phase,  are  assumed  to  vary  linearly  between  their 
c.oiturated  and  unsaturated  values  in  direct  proportion  to  the 
liquid  water  content.  If  and  K2g  are  specified  as  the 

saturated  values  of  these  coefficients,  as  given  in  Eq.  (32)  of 
Ref.  7,  then  in  the  new  translton  regime 


K 


1 


1 + 


1)5/ (h.2) 


.1/2 


(M 


(5) 


Fluctuations  in  the  saturated  value  of  humidity  are  related 
to  fluctuations  in  temperature  through  the  Clasius-Clapeynon 
equation.  For  small  fluctuations,  this  gives 


% - 


TiV: 

T R^,T 


(6) 


Thus,  in  terms  of  the  virtual  potential  temperature  which  is  used 
in  the  program 


LH 

s 

V 

R T 

m2 

V 

T 

00 

2,77? 


1/2 


1.22  + (0.61)^H' 


(7) 


This  completes  the  specification  of  the  condensation  model. 
The  influence  of  this  modification  will  be  discussed  in  connection 
v.'ith  the  sample  calculation  of  the  marine  atmospheric  boundary 
layer  in  the  next  section. 


IH 


III.  SAMPLE  CALCULATIONS 


This  section  contains  the  results  of  sample  calculations 
for  the  atmospheric  boundary  layer  over  either  land  or  water  to 
demonstrate  the  strong  difference  which  occurs  in  the  diurnal 
variations  in  these  two  types  of  boundary  layers.  The  resulting 
differences  in  the  boundary  layer  response  over  the  land  and  over 
the  water  present  the  opportunity  for  a sea-breeze  circulation 
along  any  coast.  Appendix  B presents  the  results  of  a sample 
calculation  of  this  phenomenon. 

3,1  Diurnal  Variation  Over  Land  with  Detailed  Thermal 

Radiation  Coupling 

In  previous  simulations  (2  & i})  of  the  diurnal  variation  in 
the  planetary  boundary  layer,  the  A.R.A.P,  second-order  closure 
turbulence  model  was  utilized  with  a simple  ad  hoc  radiation  loss 
term  representing  long  wavelength  cooling  of  the  boundary  layer. 
Direct  solar  heating  was  neglected.  This  radiation  loss  term  was 
assumed  known  as  a function  of  time  and  space.  We  now  present  a 
prediction  of  the  diurnal  cycle  utilizing  the  full  second-order 
closure  model  described  in  Ref.  7,  but  with  the  detailed  thermal 
radiation  model  described  in  Appendix  A.  We  consider  the  diurnal 
variation  near  relatively  dry  land  for  conditions  driven  by  the 
surface  heat  flux  measured  in  Kansas  for  a three-week  period  in 
summer  and  summarized  by  V/yngaard  (Ref.  1^1),  We  shall  compare  the 
present  calculations  utilizing  the  detailed  radiation  model  with 
the  calculation  for  the  identical  Kansas  data  reported  in  Refs,  2 
and  utilizing  an  assumed  radiation  cooling  v/hich  was  uniform  in 
time  and  represented  by  either  a uniform  or  cosine  variation  over 
space  from  the  surface  to  the  top  of  the  boundary  layer.  For  the 
present  calculation,  the  surface  humidity  mixing  ratio  is  held  at 
0,01  unless  the  saturation  value  falls  below  this  value,  in  which 
case  it  is  set  equal  to  the  local  saturated  value.  The  previous 
runs  were  independent  of  humidity. 
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The  behavior  over  the  diurnal  period  is  briefly  as  follows. 
Beginning  at  sunrise,  surface  heating  destabilizes  the  lower 
boundary  layer  and  turbulence  production  becomes  intense  over 
Che  daytime  period.  Maximum  turbulent  kinetic  energy  densities 
are  reached  at  about  2 p.m.  in  the  range  of  600-1200  m (B’lg,  3)* 
Turbulence  begins  to  subside  after  3 p.m.  in  the  afternoon  and 
has  fallen  to  a very  low  level  by  9 p.m,  in  the  evening.  Because 
of  the  viscous  coupling  with  the  upper  geostrophic  wind,  the 
coriolis-driven  wind  component  normal  to  the  geostrophic  is 
closely  coupled  to  geostrophic  during  the  day  (deviation  approxi- 
mately 6°)  and  swings  away  from  the  geostrophic  during  the 
quiescent  night  period  by  about  37°  maximum.  As  the  surface  cools 
at  night,  a low-level  inversion  in  temperature  forms  above  the 
surface.  The  cool  surface  also  induces  a slight  nocturnal  Jet. 

These  general  diurnal  phenomena  have  been  predicted  both 
by  the  model  v/ith  approximate  radiative  loss  (Ref,  2)  and  by  the 
present  model  with  detailed  thermal  radiation.  The  detailed 
radiation  model,  however,  introduces  some  quantitative  modifi- 
cations. The  most  significant  change  brought  about  by  detailed 
inclusion  of  thermal  radiation  is  in  the  diurnal  variation  in 
the  near  surface  temperature  (Pig.  ^).  With  only  an  overall 
radiative  loss  model,  the  surface  temperature  cools  bo  an  un- 
realistically low  level  during  the  night,  reaching  a maximum  of 
about  -3°C  in  potential  temperature  Just  before  sunrise.  V/ith 
detailed  inclusion  of  radiation,  the  surface  is  actually  radj- 
atively  heated  during  the  night  by  the  lower  warm,  humid  air 
layers  in  the  0-200  meter  altitudes.  The  high  humidity  of  the 
summer  midwestern  atmosphere  renders  it  optically  quite  thick  to 
thermal  radiation;  hence,  the  lower  air  layers  effectively 
blanket  the  surface  from  direct  cooling  to  space  and  lead  to  a 
much  more  realistic  minimum  in  surface  temperature  of  19°C  Just 
before  sunrise.  These  results,  of  course,  manifest  the  ff  c l 
during  the  relatively  quiescent  night  the  radiative  heat  ^ -.• 
as  important  as  the  turbulent  fluid  heat  flux. 


16 


Figure  4.  Surface  temperature  variation  as  a function  of  time  for  a typical  mldwestern 
summer  day  with  and  without  radiation  coupling. 


The  temperature  profiles  and  radiative  heabing/coollng  rate 
profiles  undergo  significant  structural  changes  over  the  diurnal 
period.  In  Pig.  5-a,  the  state  In  early  morning  before  sunrise  is 
shovm.  The  lower  level  temperature  inversion  is  near  its  maximum 
at  this  time,  and  the  boundary  layer  is  actually  being  radiatively 
heated  below  200  meters  and  radiatively  cooled  to  space  above 
200  meters.  In  Fig.  5-b,  the  state  after  sunrise  is  shown.  The 
rapid  surface  heating  has  switched  the  temperature  profile  from 
one  of  simple  inversion  to  one  of  inflection.  As  a result,  the 
lowest  surface  layers  switch  to  radiative  cooling,  the  mid-layers 
remain  in  heating,  and  the  upper  layers  remain  in  cooling.  By  mid- 
morning (Pig.  5-c),  the  inflection  is  gone  and  the  pattern  for  the 
principal  part  of  the  day  has  been  set  in  which  the  temperature 
profile  decreases  monotonlcally  from  the  surface  to  the  top  of  the 
boundary  layer,  the  lower  boundary  layer  is  in  radiative  cooling, 
and  the  upper  part  of  the  boundary  layer  is  in  weak  radiative 
heating  due  to  direct  solar  heating  slightly  exceeding  the  long 
wavelength  cooling  during  the  daytime.  In  the  early  afternoon 
with  the  sun  past  its  zenith,  the  upper  layers  switch  to  radiative 
cooling,  the  mid-layers  remain  in  heating,  and  the  lower  layers  in 
cooling  (Fig.  5-d). 

Interestingly,  the  only  period  during  the  day  in  which  the 
entire  boundary  layer  is  in  radiative  cooling  is  from  about  ^ p.m. 
in  the  afternoon  to  7 p.m.  in  the  evening.  (Recall  that  in  the 
simple  radiation  loss  model  of  previous  calculations,  the  entire 
boundary  layer  was  assumed  to  be  radiatively  cooling  at  all  points 
in  space  and  at  all  times.)  In  the  evening  (Pig.  5-e),  the  uniform 
cooling  is  about  to  give  v^ay  in  the  formation  of  the  nocturnal 
surface  inversion  which  is  well -developed  by  midnight  (Pig.  5-f) 
and  will  strengthen  to  its  maximum  by  morning  just  before  sunrise 
(Pig.  5-a). 

Detailed  predictions  of  the  radiation  field  results  in  a 
weaker  stable  temperature  profile  during  the  night.  As  a result, 
the  overall  diurnal  variation  is  somewhat  more  turbulent  with 
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Figure  5-a.  Absolute  temperature  and  radiation  flux  divergence  as  a function  of  height 
at  selected  times  during  the  day:  0600  local  time. 


3000 


Figure  5-b.  Absolute  temperature  and  radiation  flux  divergence  as  a function  of  height 
at  selected  times  during  the  day:  0900  local  time. 
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Absolute  temperature  and  radiation  flux  divergence  as  a function  of  height 
at  selected  times  during  the  day:  1A2^<  local  time. 
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Absolute  temperature  and  radiation  flux  divergence  as  a function  of  height 
at  selected  times  during  the  day:  1900  local  time. 
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Figure  5-f.  Absolute  temperature  and  radiation  flux  divergence  as  a function  of  height 
at  selected  times  during  the  day:  0200  local  time. 


detailed  radiative  prediction.  Thus,  v/ith  detailed  radiative 

prediction  the  maximum  turbulent  square^  velocity  fluctuations 

are  as  high  as  6,2  m /sec  , occur  slightly  earlier  in  the  day 

(2  p.m.),  and  at  somewhat  higher  altitude  (900  meters)  than  with 

2 2 

previous  calculations  (^,8  m /sec  maximum  at  about  4 p.m.  at 
about  600  meters). 

The  fact  that  the  nocturnal  turbulence  is  not  quite  as 
strongly  damped,  reduces  the  maximum  angle  at  which  the  surface 
wind  can  deviate  from  geostrophic  (36*^  as  opposed  to  42°)  and 
locates  the  occurrence  of  the  minimum  angle  deviation  earlier  in 
the  day  (about  10  a.m.  rather  than  noon  time).  Pig.  6.  These  re- 
sults, obtained  with  detailed  radiation  predictions,  are  closer  to 
the  wind  angle  data  of  Hoxit  (Ref.  15)  than  those  obtained  with 
the  simple  radiation  loss  term. 

3.2  Diurnal  Variation  in  the  Marine  Atmospheric  Boundary  Layer 

For  a sample  calculation  of  the  one-dimensional,  unsteady 

boundary  layer  over  the  ocean,  we  will  assume  that  the  ocean’s 

surface  temperature  is  constant  at  22°C,  that  the  geostrophic 

wind  is  constant  at  10  m/sec,  that  the  upper  level  humidity 

-3 

ratio  is  constant  at  2.5  ^ 10  , and  that  the  upper  level  po- 

tential temperature  gradient  is  0.006°C/m,  The  only  diurnal 
variation  introduced  is  in  the  solar,  short-wave  radiation, 
which  goes  through  a 12-hour  sine  wave  from  0600  local  time  to 
1800. 

Due  to  the  large  temperature  gradients  which  occur  at  cloud 

tops,  this  calculation  has  also  been  run  with  the  quasi-equilibrium 

assumption  outlined  earlier,  so  that  time  steps  are  not  limited  to 

a fraction  of  the  local  Brunt-Vaisala  period.  The  run  is  started 

with  a mixed  layer  depth  of  500  meters  and  runs  through  three  day’s 

of  simulated  time.  Although  the  results  show  a clear  diurnal 

cycle,  they  have  a substantially  larger  nonperiodic  drift  than  is 

evident  in  the  previous  results  over  land.  Figure  7 gives  the 

2 

contours  of  constant  q as  a function  of  time  and  altitude  for 


CONTOUR  q^(MVSEC^) 


TIME(HRS) 

Turbulent  kinetic  energy  as  a function  of  time  of  day  and  altitude  for  a marine 
day.  Ug  = 10  m/sec,  f = lO^^sec"*,  Tg  = 295°K,  [30/3z]  = 6®C/km.  Second  day 

of  simulation  shown.  ” 


day  two  of  the  simulation  The  growth  in  mixed  layer  depth  evi- 
dent throughout  this  day  continues  for  the  full  three-day  simu- 
lation The  diurnal  cycle  shows  up  most  clearly  in  the  lower  half 
of  the  mixed  layer,  where  the  turbulence  is  larger  during  the 
nocturnal  hours  than  during  the  afternoon  hours.  This  is  the 
opposite  condition  from  that  generally  found  over  land  (see 
Pig.  6).  The  surface  conditions  are  unstable  (negative  Monin- 
Obukhov  length)  throughout  the  day,  but  more  unstable  at  night 
due  to  a combination  of  constant  surface  temperature  and  rela- 
tively strong,  long-wavelength  cooling.  The  instability  is  re- 
duced in  the  afternoon  by  the  solar  heating  of  the  air. 

Figure  8 shows  the  contours  of  liquid  water  content.  The 
maximum  liquid  water  content  occurs  at  the  top  of  the  mixed 
layer  where  there  is  a thin  layer  of  scattered  clouds  throughout 
the  day.  The  maximum  thickness  of  clouds  occurs  near  midnight, 
but  is  a relatively  weak  maximum  compared  with  the  rest  of  the 
day  More  than  half  of  the  cloudiness  predicted  in  this  calcu- 
lation is  of  the  scattered  variety,  that  is,  the  turbulent 
fluctuations  in  H and  T permit  the  coexistence  at  the  same 
altitude  of  clouds  and  clear  sky..  Our  program  predicts  the  en- 
semble average  of  this  condition,  so  that  it  is  fair  to  think  of 
the  contours  in  Fig.  8 as  what  would  be  obtained  if  the  average 
over  a relatively  large  expanse  of  ocean  were  taken.  The 
fluctuations  in  H and  T play  a particularly  significant  role 
in  the  top  edge  of  the  cloud. 

Mixing  ratio  isolines  of  total  water  content,  vapor  plus 
liquid  are  shown  in  Fig.  9,  As  expected,  there  is  a sharp  drop- 
off in  H at  the  top  of  the  mixed  layer.  The  fact  that  water 
vapor  is  continually  fed  into  the  boundary  layer  at  the  surface 
is  the  prime  reason  for  the  continued  growth  of  the  mixed  layer. 
There  is  no  mechanism  for  v/ithdrawlng  water  from  the  boundary 
layer  in  the  current  model.  In  the  real  atmosphere  this  is  pro- 
vided generally  by  advection  to  regions  of  lower  humidity  or  by 
precipitation. 
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Plpcure  8.  Liquid  water  content  as  a function  of  time  of  day  and  altitude  for  the  conditions 
of  Figure  7. 


The  contours  of  constant  0^  are  shovm  in  Pig.  10.  The 
temperature  gradient  at  the  top  of  the  mixed  layer  is  accentuato-j 
by  the  strong  radiational  cooling  of  the  cloud  tops  which  occur 
here.  Even  when  solar  radiation  is  at  its  peak,  the  top  of  the 
clouds  are  cooled  by  long-wave  radiation.  This  is  evident  in 
Pig,  11  which  shows  contours  of  constant  radiation  flux  divergence. 
The  cloud  top  cooling  Is  sufficiently  strong  to  induce  a strongly 
unstable  temperature  gradient  in  the  upper  layer  of  the  cloud. 

This  leads  to  convectively  produced  turbulence  which  erases  the 
temperature  gradient  as  the  mixed  layer  height  is  forced  to  grow. 
This  causes  the  mixed  layer  to  grow  by  a series  of  steps  as 
evident  in  Pig.  10. 

The  wind  component  contours  in  Pigs.  12  and  13  complete  the 
description  of  the  mean  flow  variables  in  the  boundary  layer.  The 
upper  part  of  the  mixed  layer  has  slightly  supergeostrophic  values 
of  U , Values  of  V are  relatively  small  throughout  the  day. 


TIME(  HRS ) 

Figure  10.  Virtual  potential  temperature  as  a function  of  time  of  day  and  altitude  for  the 
conditions  of  Figure  7. 
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IV  CONCLUSIONS  AND  RECOMMENDATIONS 

Sample  calculations  shov;  that  the  two  major  modifications 
to  i',R.A,P.'s  marine  boundary  layer  code  during  the  past  year 
yield  reasonable  physical  results.  The  increase  in  the  di- 
mensions of  the  program  was  used  to  make  sample  calculations  of 
the  sea-breeze  circulation.  The  coupled  thermal  radiation  model 
was  used  to  make  sample  calculations  of  the  typical  diurnal 
variations  in  the  boundary  layer  over  both  land  and  water.  The 
physical  characters  of  all  of  these  results  appear  to  agree,  at 
least  qualitatively,  with  observations. 

To  increase  the  confidence  level  in  our  model,  it  will  be 
necessary  to  do  more  quantitative  verification  calculations.  Two 
ways  we  believe  this  could  be  done  over  the  coming  year  are: 

(1)  to  select  a few  well  measured  cases  from  the  GATE  (Global 
Atmospheric  Research  Program  - Atlantic  Tropical  Experiment)  data 
file  and  make  detailed  comparisons  between  our  model  predictions 
and  field  observations;  and  (2)  to  compare  our  model  results  with 
other  sophisticated  turbulent  model  calculations  such  as  have  been 
performed  in  Refs.  l6  and  17  using  a subgrid  scale  closure  scheme 
while  attempting  to  compute  all  of  the  three-dimensional,  unsteady 
turbulent  eddies  above  a given  size.  These  verification  compari- 
sons should  permit  a detailed  assessment  of  the  accuracy  of  our 
model . 

Currently,  only  the  one-dimensional  program  has  the  full 
capability  of  permitting  condensation  and  a coupled  radiation 
code.  This  capability  should  be  incorporated  into  the  two- 
dimensional  program.  This  extension  will  give  the  program  the 
capability  of  computing  coastal  fog  conditions. 

Other  modifications  to  the  program  which  it  appears  de- 
sirable to  consider  are  the  removal  of  the  hydrostatic  assumption 
from  the  two-dimensional,  unsteady  program  to  permit  features  such 
as  the  sea-breeze  gust  front  to  be  better  defined;  the  incorpo- 
ration of  the  capability  of  computing  the  flow  over  a step  change 
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in  elevation  as  occurs  at  many  shoreline  locations,  and  the 
incorporation  of  the  capability  of  predicting  the  turbulent 
diffusion  of  aerosols  of  a known  composition . 

The  ultimate  goal  of  our  program  is  to  provide  the  Navy 
with  detailed  predictions  of  the  low-level  distributions  of 
pertinent  atmospheric  variables.  V/lth  the  accomplishment  of 
the  above  tasks,  in  particular  the  detailed  verification  compari- 
sons, the  program  should  be  ready  to  test  by  investigating  distri- 
butions at  some  specific  sites  and  at  given  times  of  Interest  to 
the  Navy. 
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APPENDIX  A 


RADIANT  HEAT  TRANSFER  PREDICTION 
IN  CLEAR,  CLOUDY,  AND  TURBID 
ATMOSPHERIC  BOUNDARY  LAYERS 


D.  A.  Oliver 


I . INTRODUCTION 


A description  of  the  principal  atmospheric  radiation  effects 
in  the  planetary  boundary  layer  has  been  developed  for  incorpo- 
ration in  the  A.R.A.P.  atmospheric  boundary  layer  fluid  mechanic 
models.  This  radiation  model  includes  all  important  radiation 
effects  in  the  lov;er  atmosphere  to  at  least  a representative 
approximation.  Certain  radiative  processes,  however,  are  not 
treated  in  refined  detail  (notably,  scattering  (in  contrast  to 
absorption)  of  direct  solar  radiation  by  clouds  and  aerosols). 

The  radiation  model  is  for  the  average  one-dimensional  radiation 
field;  hence,  horizontal  gradients  and  structure  are  not  described 
except  in  an  average  sense. 

At  present,  the  statistical  coupling  of  the  radiation  and  fluid 

fields  is  not  considered.  In  principle,  the  turbulent  kinetic 

energy  destruction  v/hich  occurs  at  the  high  wave  number  end  of  the 

energy  cascade  could  be  influenced  by  radiative  transport.  Such 

effects  have  been  examined  in  stellar  atmospheres  while  very 

limited  fluid  turbulence  interaction  with  a fluctuating  radiation 

1 

field  has  been  examined  by  Goody. 

The  present  radiation  model  is  directed  to  the  interaction 
with  the  mean  fluid  equations  where  the  radiation  cooling/heating 
enters  directly  into  the  mean  energy  flow.  The  effect  on  the 
turbulence  can,  however,  be  quite  profound.  Local  radiative 
cooling  or  heating  can  alter  the  local  temperature  profile  there- 
by overturning  the  stability  of  the  fluid  or  rendering  it  stable. 
Hence,  atmospheric  radiation  can  have  an  important  coupling  with 
atmospheric  turbulence  through  the  mean  flow  and  mean  radiation 
fields.  The  present  model  is  constructed  to  reveal  this  inter- 
action. 


II.  SUMMARY  OF  THE  RADIATION  MODEL 

The  present  radiation  model  consists  of  two  principal  com- 
ponents: (1)  the  representation  of  direct  solar  heating  concen- 
trated in  the  0.7u  to  2\i  wavelength  range,  and  (2)  the 
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representation  of  terrestial  thermal  radiation  concentrated  in 
the  2\i  to  wavelength  range.  These  two  fields  are,  in 

effect,  treated  separately  in  their  interactions  with  matter 
with  respect  to  absorption,  emission,  and  scattering.  They  are 
then  combined  to  yield  the  net  heating  or  cooling  of  the  fluid. 

The  model  thus  provides  for: 

1.  Calculation  of  radiant  fluxes  and  heating  due  to  direct 
solar  radiation  including  the  effects  of: 

a.  atmospheric  molecular  absorption  by  CO2 

b.  atmospheric  molecular  absorption  by  H2O  vapor 

c.  atmospheric  absorption  by  aerosols 

d.  absorption  by  condensed  liquid  water  in  clouds  or  fog 

e.  reflectance  from  the  surface  over  which  the  boundary 
layer  flows 

2.  Calculation  of  radiant  fluxes  and  heating  due  to 
terrestial  (long  v/avelength)  thermal  radiation  including 
the  effects  of: 

a.  atmospheric  molecular  absorption  by  CO2 

b.  atmospheric  molecular  absorption  by  H2O  vapor 

c.  absorption  by  condensed  liquid  water  in  fogs  and  clouds 

3.  Parameterization  of  the  upper  atmosphere  radiation 
properties.  Since  the  model  is  primarily  intended  for 
use  in  the  lov/er  atmospheric  boundary  layer  (elevations 
of  3 km  or  lower),  the  radiation  interaction  of  the  upper 
atmosphere  with  the  boundary  layer  requires  specification. 
These  are  parameterized  in  the  present  model  by  the 
specification  of  three  average  upper  atmosphere  properties 

a.  average  cloudiness  (condensed  liquid  water  density)  of 
the  upper  atmosphere 

b.  average  humidity  of  the  upper  atmosphere  (mass  fraction 
of  v/ater  vapor) 

c.  average  turbidity  of  the  upper  atmosphere  (mass 
fraction  of  aerosols) 


III.  RADIATION  VARIABLES  AND  THE  RADIATIVE  TRANSPORT  EQUATIONS 


The  present  radiation  model  is  constructed  to  describe  the 
frequency  and  angular  averaged  radiative  intensity  and  fluxes. 
Further,  all  variables  are  assumed  to  be  a function  of  a single 
spatial  coordinate,  the  elevation  z . The  fundamental  radiation 
variables  for  either  the  direct  solar  radiation  or  the  long  v;ave- 
length  radiation  fields  are  as  f ollov/s : 


F 


Frequency  and  angular  integrated  upward  directed 
radiant  flux 


F 


q 

S 


P‘(H) 


Frequency  and  angular  integrated  downward  directed 
radiant  flux 

Radiant  heat  flux 

Radiant  source  function 

Upward  flux  at  the  base  of  the  boundary  layer  z = z^ 

Downv/ard  directed  flux  at  the  edge  of  the  upper 
atmosphere  z = H 


Downward  directed  direct  solar  radiation  at  the  edge 
of  the  upper  atmosphere 


t(z^,Z2)  Transmission  function  for  a given  radiation  field 
(either  direct  solar  or  long  v;avelength)  between 
levels  z^  and  Z2 


The  fluxes  F^(z)  , F (z)  are  the  basic  radiative  variables 
and  are  governed  by  the  upward  and  downward  transport  equations 

r 2 4 

in  integral  form  [Goody  , Ellingson  , Atwater  , Rodgers  and 
5 

Walshaw  ] 

F’’’(z)  = F'‘'(Zq)t  (z^,z)  + S(z)-  S(z^)t(z^,z) 

T (z^z)  dz^  (la) 
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Given  the  source  function  S(z)  , transmission  function  x(z^,Z2) 
and  the  boundary  values  P^(z^)  and  F~(H)  , Eqs.  (1)  uniquely 
determine  the  flux  fields  F^(z),P~(z)  . The  radiative  heat  flux 
is  then 

q^(z)  = f'^(z)  - P"(z)  (2) 

v;hile  the  radiative  heating/cooling  rate  is  given  by 

-V  • r = - - P"(z)]  (3) 

The  absorption,  scattering,  and  emissive  properties  of  the 
atmosphere  are  embedded  in  the  source  function  S(z)  and  trans- 
mission function  . These  functions  are  different  for 

the  long  wavelength  radiation  field  and  the  short  wavelength 
radiation  field.  So  also  are  the  boundary  conditions  P’^(z^) 
and  F“(H)  . 

V/e  note  that  alternative,  more  flexible  formulations  and 
solution  representations  of  the  radiative  transport  equations 
exist  [e.g.,  Marchuk”].  Hov/ever,  molecular  absorption  is  a 
principal  atmospheric  absorption  mechanism  and  the  spectra  are 
therefore  discrete  rather  than  continous.  Line  profile  dis- 
tortion (variation  of  absorption  over  the  line)  has  a profound 
effect  on  atmospheric  radiation  absorption  and  cannot  be  readily 
described  by  schemes  such  as  those  of  Marchuk^.  Neglect  of 
line  profile  distortion  for  molecular  absorbers  like  CO^  and 
water  vapor  can  lead  to  errors  of  the  order  of  30%  in  absorption 
prediction.  The  integral  formulation  given  above  (though  less 
flexible  with  re'spect  to  boundary  conditions)  can  readily  in- 
corporate the  line  profile  distortion  characteristic  of  molecular 
absorbers  in  the  transmission  function. 


IV.  TERRESTIAL  RADIATION 


•f  — 

We  now  discuss  the  boundary  values  (P  (z^),F  (H))  , source 
function  (S(z))  , and  transmission  function  t(z^jZ2)  appropriate 
for  terrestial  radiation  in  the  lower  atmosphere.  Although  there 
is  evidence  that  strong  nonequilibrium  effects  exist  in  the  upper 
rarefied  atmosphere,  Kirchhoff's  lav;  holds  quite  well  in  the 
boundary  layer  region  and  it  is  permissible  to  set  the  radiant 
source  function  equal  to  the  Planck  function: 

S(z)  = B(z)  = ot\z)  (^) 

(a  is  the  Stephan-Boltzmann  Constant). 

Appropriate  vlaues  for  the  boundary  conditions  are 

P^Zo)  = B^  (5) 

P"(H)  = eB(H)  (6) 

In  the  above,  B^  is  a boundary  Planck  function  which  in  general 
is  equal  to  B(z^)  ; however,  a discontinuity  between  surface 
and  lower  atmospheric  temperature  is  permitted  for  certain  kinds 
of  modelling  in  which  the  analyst  need  not  specify  B^  = B(z^)  . 

The  parameter  e is  the  upper  atmosphere  emissivity  for 
downward  radiation.  If  z = H lies  outside  the  atmosphere  and 
no  thermal  radiation  is  incident  upon  the  earth,  a rigorous  con- 
dition is  e = 0 . In  practice,  it  is  found  that  a value  of 
e ~ 0.1  - 0.15  applied  at  the  "edge"  of  the  upper  atmosphere 
(p  ~ 10  mb)  gives  reasonable  agreement  with  experimental  measure- 
ments of  F~(z)  [Ellingson'^]. 

The  transmission  function  xCz^jZg)  describes  the  atmospheric 
interaction  with  the  radiation  field.  For  the  three  components 
which  interact,  let  us  denote  by  subscripts  i = 1,2,3  the 
following. 

p^(z)  = mass  density  of  H2O  vapor 
P2(z)  = mass  density  of  CO2 
P^{z)  = mass  density  of  H2O  liquid 
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iii.iiWin^ir 


The  transmission  function  x (2^,22)  can  then  be  considered  to 
be  a function  of  the  mass  areal  densities  m^^  in  the  interval 
(2^,22)  as 

t(22^322)  *“  x(m^^m2>m^)  *' 


v/here  the  mass  areal  density  for  species  i , m^^  is  defined  as 


The  transmission  function  for  long  wavelength  radiation  is  modelled 

O 

after  that  of  Feigel'son  : 

r ^ 

T(m^,m2,m^)  = T2(m2)|  Y1  (-Pot^^^m^) 

^k=l 

+ TT  Y1 

n=2  k=l  K n 1 ^ 

The  function  T^(m2)  accounts  for  cloudiness: 

T^(m2)  = exp  [-Bam^]  (9) 

In  the  above  transmission  functions  the  coefficient  B is  the 

diffusivity  factor  which  accounts  for  the  fact  that  the  terrestial 

g 

radiation  is  diffu^’e  rather  than  collimated  [B  = 1.66  , Elsasser  ]. 
The  coefficients  ^k^^  > ^k”^  have  been  computed  by  Peigel’son^^ 
[Table  III. 2. 3,  p.  83]  to  accurately  match  the  summary  of  CO2  and 

H^O  vapor  laboratory  measurement  data  compiled  by  Davis  and  Vie2ee 

‘^2  2 
For  m^  < 10  gm/cm  and  m2  < 100  gm/cm  the  Feigel’son  trans- 
mission function  matches  the  data  summary  of  Davis  and  Viesee  within 
5%  for  CO2  and  H2O  vapor. 

The  coefficient  a in  the  transmission  component  T^  accounts 
for  the  absorptive  properties  of  liquid  water.  For  long  wave- 
length thermal  radiation  in  the  range  from  2 to  ^Ou  the  absorptive 
characteristics  of  water  droplets  depend  strongly  upon  their  sise. 
Hence,  an  accurate  description  of  the  absorption  in  a cloud 


I 


r; 
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requires  specification  of  the  droplet  size  distribution  within 
the  cloud. 

12  11 

Feigel'son  and  Gradus,  et  al.  have  determined  the 
effective  absorption  coefficient  a for  droplet  radii  distri- 
bution functions  of  the  form 


N(r)  = N^r*^e“^^ 

where  (n,a)  are  parameters  characterizing  the  distribution. 

These  parameters  are  in  turn  related  to  cloud  type.  Small  values 
of  n and  a correspond  to  so-called  "narrow"  distributions, 
while  large  (n,a)  correspond  to  "wide"  distributions.  Narrow 
distributions  are  characteristic  of  Ns  , St  , II  , while  wide 
distributions  are  characteristic  of  Sc  , As  . Clearly,  wide 
distributions  are  clouds  dominated  by  large  droplets  while  narrow 
distributions  are  dominated  by  small  droplets.  Fogs  or  low-lying 
marine  clouds  would  perhaps  be  more  typical  of  the  narrow  distri- 
bution. These-  distribution  functions  are  utilized  with  the 
spectral  absorption  cross  sections  of  droplets  of  different  radii 

ill  1 15 

tabulated  by  Herman  and  Zel'manovich,  et  al.  to  yield  the 
total  absorption  coefficient  a . The  results  of  these  detailed 
calculations  for 

narrow  distribution 


wide  distribution 


and 


n 

a 

■ ave 

n 

a 

*1 

ave 


are 


C 2 

I 500  cm  /gm  (wide) 

I 1700  cm^/gm  (narrow) 
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The  factor  of  3 or  more  difference  between  these  two  absorption 
coefficients  points  out  clearly  that  accurate  prediction  of 
radiation  in  a cloudy  atmosphere  requires  a correlation  between 
cloud  type  (n,a,r^^g)  and  its  local  thermodynamic  state 
Temperature),  and  perhaps  its  condensation  history.  Such  detail 
is  not  included  in  the  present  radiation  model.  Rather  we  treat 

the  absorption  a as  a cloud  type  parameter  varying  between  500 
2 

and  1700  cm  /gm. 


V.  DIRECT  SOLAR  RADIATION 


+ — 

We  nov;  turn  to  the  boundary  values  [P  (z^),F  (H)],  source 
function  S(z)  and  transmission  function  xCz^jZg)  appropriate 
for  direct  solar  radiation.  The  approximation  that  emission  into 
the  short  wavelength  spectrum  is  negligible  is  made  so  that  for 
the  direct  solar  radiation  field  we  set 

S(z)  = 0 

The  boundary  values  are  set  as 

^■(H)  = Ig 

= rp'(z^) 

where  f is  the  surface  reflectance  for  direct  solar  radiation. 
The  short  wavelength  transmission  function  formally  depends  upon 
four  absorbing  areal  mass  densities  as 

t(z^,Z2>  = T(m^,m2,m2,ni^) 

where  m^^Ci  = 1,2,3)  are  as  defined  in  Eq.  (7)  and  m^^  is  the 
aerosol  areal  mass  density  defined  in  terms  of  the  volum'^trlc 
aerosol  density  as 


zi 


The  short  wavelength  transmission  function  is  then  given  by 

T(m^,m2,m2,m^|)  = T^(m^)  • • T^Cm^)  • (H) 

Band  overlap  between  CO2  and  H2O  vapor  can  be  neglected  for 
the  short  wavelength  radiation;  hence,  the  total  transmission 
function  can  be  expressed  as  the  product  of  the  individual  specie 
transmission  functions.  Further,  line  profile  distortion  is 
important  for  the  water  vapor  transmission  function  T^  so  we 
represent  it  in  the  functional  form 


(12) 


T^(tn^)  = 1 


) 


which  embraces  the  functional  forms  for  strong  and  weak  band 
absorption  [Sivkov^^].  The  parameters  a^  , b^  are  selected 
as 

a = 0.04 

V 

\ = 15 
2 

when  m^  is  in  units  of  gm/cm  . The  water  vapor  transmission 

function,  Eq.  (12),  with  these  parameter  choices  of  a^  j by 

matches  the  data  of  Sivkovl"^  and  Manabe  and  Strickler^®  for 
2 

m^  < 3 gm/cm  within  5/5  • 


The  transmission  functions  for  v'ater  vapor  free  air  and 
liquid  water  are  represented  as 


Tg(m2)  = exp  [-Sa^m^]  (13) 

T^(m^)  = exp  [-(jajj^m^]  (l4) 


The  absorption  coefficient  for  the  gaseous  molecular  components 
of  clean  water-free  air  is  taken  as 


-4  2 

a = 0.7  X 10  cm  /gm 
6 

which  matches  well  the  data  of  Sivkovl*^,  Manabe  and  Stricklerl®, 

iq 

and  Eltermann,  et  al. 

In  the  case  of  liquid  water  absorption  of  direct  solar 

20 

radiation,  Peigel'son  notes  that  for  wavelengths  less  than  2y 

the  absorption  is  not  sensitive  to  droplet  size.  Further,  the 

liquid  water  absorption  spectrum  is  continuous  so  that  the  form 

Eq.  (13)  is  quite  accurate.  Integrating  the  spectral  liquid  water 

21 

absorption  coefficient  tabulated  by  Feigel’son  over  the  direct 
solar  spectrum  we  find 

2 

= 13  cm  /gm 
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Lastly,  v;e  consider  the  transmission  Tg^(m|j)  of  the  aerosol. 
Since  the  absorotion  by  qerosol  particles,  like  that  of  water 
droplets,  can  be  considered  to  be  continuous,  the  transmission 
function  may  be  represented  as 

Ta(m,j)  = exp  C-a^m^]  (15) 

The  absorption  coefficient  a„  is  a sensitive  function  of  the 

“ 22 

particle  size  distribution.  As  noted  by  Kondratyev  , the  particle 
size  distribution  of  typical  atmospheric  aerosols  varies  signifi- 
cantly with  height;  hence,  a detailed  description  of  radiative 
interaction  with  the  aerosol  requires  specification  (or  simultaneous 
prediction)  of  the  particle  size  distribution  function  as  well  as 
the  total  aerosol  mass  density.  In  the  event  that  and 

are  predictable  or  known  in  such  detail,  the  form  Eq.  (15)  is 
appropriate  and  useful. 

In  lieu  of  a detailed  prediction  procedure  for  and  , 

three  approximate  methods  varying  in  their  degree  of  detail  are 
now  presented  for  the  representation  of  the  aerosol  absorption. 

2 •2  2k 

For  (globally)  average  conditions  Ivlev  ^ and  Toon  and  Pollock 
have  summarized  global  aerosol  concentrations,  particle  distri- 
butions, and  absorption  coefficients  (or  optical  depths  at  0.55P 
Tq  ^^).  For  the  boundary  layer  region  (0  < z < 5km)  the  Ivlev 
and  Toon-Pollock  data  are  summarized  in  Table  I.  The  absorption 
data  of  Toon  and  Pollock  is  systematically  below  that  of  Ivlev  at 
the  lower  altitudes.  Perhaps  this  is  in  part  because  it  is  not 
averaged  over  frequency  but  centered  at  0.55u  . 

The  first  and  simplest  average  treatment  of  aerosols  is 
based  on  the  assumption  that  local  absorption  is  equal  to  the 
global  average.  The  transmission  function  T.  is  then  repre- 

CL 

sented  as 

T^  = exp  C-a*z] 

v/lth  a*  taken  from  Table  I. 
a 


(16) 


nmm 


I 

I 

I 

I 


TABLE  I.  SUMMARY  OP  PARTICLE  CONCENTRATION  AND 
ABSORPTION  DATA  OF  IVLEV^^  AND  TOON  AND  POLLOCK^^ 


; 

z(km) 

Particle  Concentration 

aj(km"^) 

a*(km“^) 

(Ivlev) 

(cm“^) 

(Prom  Toon  and 
Pollock  calc, 
from  Tq  55) 

(Ivlev) 

r < 0.1 

2r  > 0.1 

0-1 

h 

2x10^ 

1.5x10^ 

il,2xl0“^ 

1.2x10"^ 

1-2 

3><10^ 

8.0x10^ 

2.5x10“^ 

6x10"^ 

j*  >- 

2-3 

IxlO^ 

2.7x10^ 

1.6x10“^ 

3x10"^ 

* 

5x10^ 

2.0x10^ 

0.8x10"^ 

1.8x10"^ 

I|-5 

3x10^ 

1.5x10^ 

6x10"^ 

7x10"^ 

'f 
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A second  approach  which  allows  some  detail  in  the  variation 
of  aerosol  concentration  is  the  following.  The  transmission 
function  is  represented  in  terms  of  the  particle  density  N(z) 
rather  than  the  mass  density  P|j  as 

Tg^(na)  = exp 


where 


n. 


2 

N( z) dz 


(18) 


From  the  data  of  Table  I,  the  absorption  coefficient  a may  be 

cL 

computed  and  is  tabulated  in  Table  II. 


TABLE  II.  SPECIFIC  PARTICLE  ABSORPTION 

COEFFICIENT  a COMPUTED  FROM  THE  DATA  OF  TABLE  I 
a 


z(km) 

^ / 2 \ 
a^iom  ) 

0-1 

6x10"^ 

1-2 

5.2x10“^ 

2-3 

3.75x10"^ 

2.9*<xl0"^ 

i»-5 

l.i<lxl0“^ 

The  representation  in  terms  of  would  be  useful  in  situations 

where  the  particle  size  distribution  could  be  assumed  invariant 
with  respect  to  total  particle  density  N . One  would  expect 
this  to  be  the  case  for  low  aerosol  densities  where  aerosol 
particle-particle  interactions  are  rare.  Aerosol  absorption  would 
then  be  a function  only  of  aerosol  concentration  N(z)  with  the 
absorption  coefficient  taken  to  possess  a universal  height 

O. 

distribution  typical  of  the  global  average  of  Table  II.  The 


number  density  N(z)  Is  then  predicted  with  a turbulent  species 
transport  equation  as  is  done  for  water  vapor  or  pollutants  in 
existing  A.R.A.P.  models.  Note  also  that  the  absorption  co-  . 

efficient  a„  could  be,  within  the  limits  of  the  approximations 

“ ~ 2 
Involved,  taken  as  a constant  = 5 10  cm  in  the  0-3  km 

ci 

region  which  is  the  boundary  layer  region  of  principal  interest. 

A third  approach  similar  to  the  second  in  level  of  detail 
represents  the  transmission  function  as 

t(Zj^.Zj)  = T^(m^)  • ymp  • T^(r,3) 

where  m^  , m^  are  as  defined  in  Eq.  (7),  but  m^  is  constructed 
to  be  a composite  of  clear  air  molecular  absorption  and  aerosol 
absorption  in  the  following  manner: 

^2  ~J  (19) 

^1 

In  Eq.  (19),  p is  the  air  density,  n(z)  is  the  mass  fraction 

s 

of  aerosol  relative  to  the  mass  fraction  at  some  reference 
elevation  (typically  the  surface).  The  parameter  t^  is  the 
reference  turbidity  defined  as  the  ratio  of  the  absorption  co- 
efficient of  the  aerosol  to  the  absorption  coefficient  of  clear 
air  at  the  reference  level.  Thus,  in  this  approach  the  aerosol 
transmission  is  represented  in  terms  of  two  parameters:  a non- 
dimensional  turbidity  profile  n(z)  and  a reference  turbidity 
level  t^  . Reasonable  assumptions  may  be  made  for  these  para- 
meters for  specific  calculational  cases.  For  example,  in  the 
marine  boundary  layer,  the  turbidity  primarily  arises  from  sea 
salt  particles  turbulently  diffused  and  convected  from  the  surface 
as  is  the  water  vapor.  In  this  case,  one  may  make  the  turbidity 
profile  n(z)  similar  to  the  humidity  profile. 


25 

For  values  of  reference  turbidity,  Sivkov  proposes  the 
following  semi-quantitative  classification: 

t^  Atmospheric  Condition 

0 free  of  aerosols 

1 mildly  turbid  atmosphere 

4 highly  turbid  atmosphere 
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VI.  REPRESENTATION  OF  THE  UPPER  ATMOSPHERE 


The  mass  integrals  m^^  defined  by  Eq,.  (7)  extend  through- 
out the  entire  atmosphere.  For  the  domain  of  the  mass  integral 
v/hich  extends  above  the  boundary  layer  height  z = H^^  some 
assumption  about  the  upper  atmosphere  properties  are  required. 
These  upper  atmosphere  characteristics  are  represented  by  three 
model  parameters: 


t^  = average  turbidity  of  the  upper  atmosphere 

h^  = average  humidity  of  the  upper  atmosphere 

= average  liquid  water  fraction  of  the  upper  atmosphere 
(If  the  upper  atmosphere  is  cloudless,  = 0 , 
except  for  ice  clouds) 


Let  H^  be  the  atmosphere  scale  height  (H^  = 799^m) . Then  for 
the  mass  integrals  that  extend  into  the  upper  atmosphere,  we  have 


Hb 

J * (">uh 


where  is  defined  as 


and 


p^(z^)AH 


AH  = H^-  H^^ 


In  particular,  for  the  four  mass  integrals  we  have 


(m^)2 


<^>2 


Pj^(Zo)h„AH 

p^(z^)AH 


VII.  MODEL  COMPARISON  WITH  ATMOSPHERIC  MEASUREMENTS 


For  the  purpose  of  initially  verifying  the  general  correct- 
ness of  the  present  radiation  model,  we  have  selected  a comparison 

of  long  wavelength  radiation  data  reported  from  the  International 

2 S 

Radiometersonde  Intercomparison  Program  by  Gille  and  Kuhn‘S  and 

a verification  of  short  wavelength  solar  heating  predictions  of 

the  model  compared  to  the  observations  of  Reynolds,  Yonder  Haar, 

27 

and  Cox  ' in  the  BOMEX  marine  boundary  layer. 

In  the  case  of  the  Gille-Kuhn  Intercomparison  data,  at- 
mospheric temperature,  pressure,  and  humidity  as  measured  are  used 
as  inputs  to  the  model  from  which  the  long  wavelength  radiant 
fields  and  cooling  rates  are  predicted.  These  are  then  compared 
with  the  measured  radiant  fields.  The  results  of  this  comparison 
for  two  typical  cases  are  shown  in  Pigs.  1-6,  one  for  clear  skies 
and  one  with  high  overcast  in  the  upper  atmosphere.  It  can  be 
seen  that  the  comparison  between  predicted  and  experimental  values 
is  reasonable. 

As  a check  on  the  direct  solar  heating  predictions  of  the  model, 

we  have  computed  the  direct  solar  heating  rate  in  the  lower  marine 

boundary  layer  for  conditions  generally  similar  (but  not  identical 

in  detail)  to  that  in  the  BOMEX  experiment  reported  by  Reynolds, 

27 

Yonder  Haar,  and  Cox  . The  boundary  layer  was  assumed  to  have  a 
reference  turbidity  of  t^  = 1 (mildly  turbid)  and  the  heating- 
rates  were  averaged  over  daylight  hours  to  be  consistent  with  the 
normalization  utilized  by  Reynolds,  et  al.  The  model  predicted 
heating  rate  is  compared  vjith  the  measurements  in  Fig.  7.  A 
favorable  comparison  of  order  of  magnitude  again  results. 
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VIII.  AN  APPROXIMATE  FORM  OF  THE  TRANSFER  EQUATIONS 


The  transfer  Eqs.  (1)  contain  a convolution  of  the  source 
function  gradient  over  the  transmission  function  of  the  form 

^2 

For  practical  calculations  the  computation  of  this  convolution 
Integral  can  be  extremely  lengthy.  An  approximate  but  rapid 
computation  would  therefore  be  desirable.  The  transmission 
function  t(z",z)  Is  a monotone  decreasing  function  of  the 
difference  of  arguments  z - vanishes  for  z - z^  , 
and  Is  Identically  unity  for  z'  = z . It  Is  easy  to  show 
that  there  exists  a suitably  defined  average  < > such  that 


<TVZ 


[SCZj)  - 


SCZ]^)  ] 


(20) 


As  an  approximation  to  this  average,  we  shall  adopt  a simple 
definition  of  the  average  transmission  2 as  the  trans- 

mission function  calculated  In  terms  of  the  arithmetic  average  of 
the  argument  difference,  viz: 


<t(z%z)>„  „ = T(m  /2, 

^1^  "2 


m2/2. 


m^/2) 


(21) 


where  the  m^  are  defined  as  In  Eq.  (7).  The  approximation, 

Eq.  (21),  of  the  convolution  Is  computationally  extremely  rapid. 
A comparison  of  the  approximate  form,  Eq.  (20),  with  the  exact 
convolution  form,  Eqs,  (1),  Is  shown  In  Figs.  8-10. 
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IX.  CONCLUSION 


We  conclude  that  the  radiation  model  we  have  developed  for 
incorporation  in  A.R.A.P.  atmospheric  fluid  dynamic  calculations 
yields  results  which  are  consistent  with  the  limited  experimental 
radiation  measurements  available.  Further  confirmation  of  its 
validity  will,  of  course,  come  in  the  use  of  the  full  fluid- 
radiative  atmospheric  equations  in  predicting  marine  boundary 
layer  phenomena  of  general  interest. 
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Figure  A-1.  Upvmrd  long  wavelength  radiant  flux  in  the  earth's 
atmosphere  corresponding  to  Case  D-3,  clear  skies, 
J,C.  Gille  and  P.M.  Kuhn^o  compared  to  radiation 
model  predictions. 
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Figure  A- 


. Comparison  of  radiation  model  predictions  and  experimental 
summary  of  J.C.  Gille  and  P.M.  Kuhn^°  of  atmospheric  long 
v/avelength  cooling  rate.  Case  D-3,  clear  skies. 
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Figure 


-5.  Downward  long  wavelength  radiant  flux.  Case  D-5  (overcast) 
of  J.C.  Gille  and  P.M.  Kuhn^®  compared  to  radiation 
model  predictions. 
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A-6.  Comparison  of  radiation  model  predictions  and  experimental 
summary  of  J.C.  Gille  and  P.M.  Kuhn^o  of  atmospheric  long 
wavelength  cooling  rate.  Case  D-5  (overcast). 
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Figure  A-7.  Comparison  of  direct  solar  heating  rate  predicted  by 
the  radiation  model  and  measured  heating  rates  of 
Reynolds,  et  al.^'.  Bulge  in  the  experimental  data 
is  due  to  wind-borne  Sahara  dust  between  2000-4000  m. 
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Comparison  of  approximate  and  exact  downward  radiant  | 

flux  for  conditions  of  Fig.  8.  j 
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Figure 


-10.  Comparison  of  approximate  and  exact  cooling  rates  for 
conditions  of  Fig.  8. 
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A SECOND-ORDER  CLOSURE  MODEL  OP  TURBULENT  TRANSPORT 
IN  THE  COASTAL  PLANETARY  BOUNDARY  LAYER* 


W.  S.  Lewellen  and  M.  E.  Teske 

Aeronautical  Research  Associates  of  Princeton,  Inc. 
Princeton,  New  Jersey 


INTRODUCTION 

The  diurnal  variation  In  the  plane- 
tary boundary  layer  Is  quite  different 
over  land  than  It  is  over  a large  body  of 
water.  This  variation  leads  to  the 
familiar  sea-breeze  circulation  pattern 
near  the  shore.  Many  of  the  character- 
istics of  lake/sea  breezes  have  been  de- 
scribed In  a review  by  Lyons  (1975)  with 
particular  emphasis  on  pollutant  trans- 
port In  such  an  environment.  The  most 
recent  attempts  to  model  this  phenomenon 
have  been  the  work  of  Plelke  and  Mahrer 
(1975)  and  Mak  and  Walsh  (1976).  In 
order  to  simulate  the  stability  variation 
of  the  turbulent  transport  processes,  an 
essential  part  of  the  sea-breeze  phenome- 
non, both  works  use  an  empirically 
assumed  eddy  viscosity  that  Is  a function 
of  both  space  and  time. 

Turbulent  models  based  on  second- 
order  closur  ”s  proving  to  be  very  use- 
ful for  solvi  ^ turbulent  transport  of 
momentum,  energy,  and  species  In  the 
planetary  boundary  layer  (Lewellen  and 
Teske,  1975;  Wyngaard,  et  al.,  197*1; 

Mellor  and  Yamada,  197*1).  The  purpose  of 
this  paper  Is  to  present  some  prelininary 
results  of  a calculation  of  the  sea-breeze 
circulation  using  a second-order  closure 
model  of  turbulence.  Such  an  approach 
eliminates  the  need  to  specify  an  empiri- 
cal eddy  viscosity  variation. 

REVIEW  OF  SECOND-ORDER  CLOSURE  MODEL 

Work  at  A.R.A.P.  over  the  last  few 
years  has  led  to  the  development  of  an 
Invariant,  scuond-order  closure  technique 
for  the  determination  of  the  Reynolds 
stress  terms  Uj[U^  , and  the  turbulent  flux 

of  potential  temperature u^e  and  flur.  of 

humidity  mixing  ratio  mK  ,‘that  appear  in 
the  mean  momentum,  energy,  and  humidity 
equations.  By  manipulating  the  fluctu- 
ating equations  of  motion,  one  can  derive 
a set  of  exact  equations  for  the  needed 
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second-order  correlations  and  the  further 
correlations  P"  , FIT  and  that 
appear  (Donaldson,  1973).  The  task  of 
second-order  closure  Is  to  model  the  re- 
sulting terms  Involving  third-order  corre- 
lations and  other  unknown  second-order 
correlations  as  functions  of  the  known 
second-order  correlations,  mean  flow 
gradients,  and  proper  scaling  and  m.odellng 
constants,  and  thus  close  the  equation  set. 
To  do  this  we  choose  the  simplest  models 
that  have  proper  tensor  symmetry,  dl- 
menslonalizatlon,  and  desired  physical 
properties;  make  a consistent  selection  of 
the  values  of  the  modeling  constants;  then 
leave  these  values  fixed  when  computing  new 
flow  geometries.  Lewellen  and  Teske  (1975) 
contains  an  extensive  discussion  of  the 
modeling  procedure  and  an  examination  of 
the  many  flow  geometries  used  to  evaluate 
and  verify  the  selection  of  modeling  con- 
stants. Modeling  predictions  for  more 
complicated  flows  Include  the  diurnal 
planetary  boundary  layer  (Lewellen,  et  al., 
197*1),  stratified  flows  (Lewellen,  et  al., 
1976).,  and  vortex  wakes  (Bllanln,  et  al., 
1976). 

A solution  of  the  full  set  of  equations 
(21  Including  a differential  equation  for 
the  turbulent  scale  length)  Is  a formidable 
numerical  problem.  For  the  present  calcu- 
lation, we  elect  to  make  a simplifying 
assumption  to  the  turbulent  correlation 
equations  by  treating  them  In  a quasi- 
steady,  nondiffuslve  limit.  The  partial 
differential  equations  for  the  second- 
order  correlations  In  this  boundary  layer 
problem  then  reduce  to  algebraic  ex- 
pressions as  a function  of  the  vertical 
gradients  of  the  two  mean  horizontal  wind 
components,  the  vertical  potential  tempera- 
ture, and  the  humidity  mixing  ratio.  The 
full  dynamic  behavior  Is  also  retained  for 
the  turbulent  kinetic  energy  q*  ap.d  the 
turbulent  scale  A . The  resulting  alge- 
braic equations  for  the  1*1  turbulent 
correlations  may  be  solved  explicitly 
(Lewellen,  et  al.,  197*0. 

The  calculations  are  performed  using 
a computer  code  that  solves  the  two-di- 
mensional, unsteady  equations  of  motion 
for  the  moan  velocities  U and  V , 
virtual  temperature  Oy  , humidity  H , 
turbulent  kinetic  energy  q*  and  turbulent 


Bcale  length  A . Turbulent  correlatlona 
uiuj,  uj^Oy,  0^  , uXIT  , Oyti  and  Fn"  are 
found  by  applying  the  quasl-equlllbrium 
approximation  discussed  earlier.  Details 
of  the  finite  differencing  technique  em- 
ployed and  the  Implicit  solution  techni- 
que used  may  be  found  In  Lewellen  and 
Teske  (1975). 

Since  the  vertical  velocity  is  of 
higher  order  In  the  planetary  boundary 
layer,  the  vertical  momentum  equation  re- 
duces to  the  hydrostatic  approximation. 

The  vertical  velocity  is  then  determined 
from  continuity  via  a stream  function- 
vortlclty  approach,  with  the  vortlclty 
determined  directly  from  the  horizontal 
momentum  equation.  The  Poisson  equation 
for  i|i  Is  solved  directly  by  an  appli- 
cation of  the  DLKTRI  general  elliptic 
solver  developed  by  Swarztrauber  and 
Sweet  (1975). 

The  advectlon  terns  In  the  present 
calculation  are  upwind  differenced.  This 
undoubtedly  adds  horizontal  diffusion  to 
the  calculation.  However,  the  principal 
balance  In  the  present  calculation  Is  be- 
tween vertical  diffusion  and  horizontal 
advectlon  which  should  dominate  the  hori- 
zontal numerical  diffusion.  The  character 
of  the  present  result  should  be  essential- 
ly correct.  The  largest  errors  Intro- 
duced are  probably  due  to  the  relatively 
large  Ay  grid  spacing  which  varies  from 
3 to  18  km, 

FORMULATION  OF  THE  SAMPLE  CALCULATION 


air,  we  have  also  Included  a thermal 
radiation  source  term  in  the  energy 
equation  equal  tc 

Qp  - -2  »«  10"'*(®C/sec)  exp(-z/500  m)  (2) 

In  a real  boundary  layer  this  radiation 
flux  divergence  term  Is  strongly  coupled 
to  the  flow  variables,  particularly  the 
water  vapor  distribution.  In  this  pre- 
liminary model,  the  simple  cooling 
function  In  Eq.  (2)  Is  Incorporated  to 
permit  the  boundary  layer  over  the  water 
surface  to  maintain  a slight  Instability 
as  observed  over  the  ocean. 

The  computation  Is  performed  over  a 
domain  stretching  from  126  km  on  either 
side  of  the  shoreline,  with  a vertical 
celling  at  1700  m.  One  of  the  principal 
difficulties  of  the  computation  Is  the 
appropriate  treatment  of  the  internal 
gravity  waves  at  the  open  boundaries.  In 
this  calculation  horizontal  diffusion  is 
Increased  smoothly  as  the  ' oundary  Is 
approached  to  damp  out  any  waves  and  pre- 
vent reflection  back  Into  the  domain  of 
Interest.  The  results  shown  below  are 
confined  to  80  km  on  either  side  of  the 
shorell  le,  since  the  wave  liner  influences 
the  results  between  80  and  126  km. 

The  geostrophlc  wind  Is  aligned  with 
the  coastline  and  has  a value  of  either 
*5  m/see.  With  the  land  surface  on  the 
right  In  the  figures  shown,  the  plus  sign 
corresponds  to  flow  out  of  the  paper. 
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In  the  present  calculation  the  land 
and  water  surface  temperature  is  taken 
equal  to  290«K  at  t “ 0 . The  land  sur- 
face temperature  Is  then  permitted  to  add 
a perturbation  which  varies  with  time  as 

«LS  ■ 

while  the  water  surface  temperature  is 
held  fixed.  Equation  (1)  Ignores  any 
spatial  variation  of  the  land  surface 
temperature.  To  account  for  this  vari- 
ation correctly,  it  would  be  necessary  to 
solve  for  heat  flow  In  the  top  layer  of 
the  land  soil  in  a completely  coupled 
manner  with  the  boundary  layer  calcu- 
lation. 

The  surface  temperature  of  shallow 
coastal  water  may  also  have  a slight  di- 
urnal variation.  However,  we  believe 
this  ideal  problem  of  a land  mass  under- 
going a sine  v/ave  variation  in  Its  surface 
temperature,  bounded  by  a constant 
temperature  v;ater  reservoir,  exhibits  the 
essential  features  of  the  diurnal  vari- 
ation In  the  sea-breeze  circulation.  We 
have  Included  a surface  roughness  dif- 
ference between  the  land  and  the  water 
with  Zq  » 0.1m  for  the  land  and 
2.3  10"’  m for  the  water. 

In  partial  simulation  of  Infrared 
radiation  cooling  from  the  moist  marine 
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The  other  variables  which  signifi- 
cantly affect  the  results  are  the  Initial  \ 

height  of  the  Inversion  layer  and  the  : 

lapse  rate  above  the  boundary  layer.  ! 

Above  z « 1100  m the  vertical  potential  j 

temperature  gradient  Is  set  at  -0.001®C/m  i 

Initially.  j 

RESULTS  ; 

\ 

Figure  1 shows  the  time  history  of  ! 

the  two  components  of  horizontal  velocity 
at  the  shoreline  (y  * 0)  as  a function  of 
altitude.  The  component  perpendicular  to 
the  shoreline,  V , Is  directed  toward  the 
land  (positive)  below  approximately  ^400  m,  ; 

except  for  a brief  period  near  the  end  of  » 

the  cycle.  The  return  flow  at  altitude  | 

is  substantially  reduced  in  magnitude  but 
occurs  over  a larger  region.  The  sea 
breeze  velocity  Increases  as  the  land  sur- 
face temperature  Increases,  with  its  maxi- 
mum value  occurring  about  an  hour  after  I 

the  maximum  surface  temperature  Is  reached  I 

at  6 hours.  The  reduction  in  the  sea  I 

breeze  lags  far  behind  the  decrease  in  the  f 

surface  temperature  so  that  very  little  I 

reversal  or  land  breeze  occurs  (only  near  \ 

21  hours).  The  velocity  component  s' 

parallel  to  the  shoreline,  U , roaches  a ] 

maximum  value  near  the  point  (12  hours)  > 1 

where  the  land  surface  temperature  goes  s | 

through  Its  zero  point.  | | 
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Pig.  1.  Time  histories  of  the  horizontal  wind  components  at  the  shoreline  (y  • 0)  for 
Uff  * 5 m/sec:  (a)  V velocity  normal  to  the  shore;  (b)  U velocity  parallel  to  the 
snore.  A positive  value  of  V Is  sea  breeze. 
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The  pressure  gradient  Induced  by  the 
temperature  difference  In  the  boundary 
layer  over  the  land  and  that  over  the 
water  appears  to  first  yield  a wind  ve- 
locity down  the  gradient,  and  then  the 
wind  vector  swings  around  to  be  perpen- 
dicular to  the  pressure  gradient  in  the 
same  manner  as  the  geostrophlc  wind 
balancing  the  geostrophlc  pressure  gradi- 
ent. This  gradual  turning  of  the  di- 
rection of  the  wind  in  the  course  of  the 
diurnal  cycle  Is  present  In  observations 
(Lyons,  1975)  and  appears  In  previous 
models  (e.g.,  Neumann  and  Mahrer,  1971). 

When  the  geostrophlc  wind  Is  re- 
versed, with  other  conditions  unchanged, 
the  maximum  value  of  the  sea  breeze  is  re- 
duced substantially,  and  a land  breeze 
occurs  for  approximately  one-fourth  of  the 
cycle  as  seen  In  Fig.  2.  The  parallel  ve- 
locity component  completely  changes  Its 
character  and  even  shows  a slight  flow  re- 
versal near  15  hours*  Although  the  geo- 
strophlc wind  Is  aligned  with  the  shore-  . 
line  In  both  cases,  the  Ekman-type  spiral 
of  the  wind  vector  In  the  layer  tends  to 
add  to  the  sea  breeze  in  Fig.  1 and  sub- 
tract (or  add  to  the  land  breeze)  In 
Fig.  2.  The  times  and  altitudes  at  which 
the  wind  parallel  to  the  shoreline  Is 
amplified  or  retarded  are  Interchanged  be- 
tween Figs,  lb  and  2b. 

The  spatial  variation  of  the  sea 
breeze  circulation  at  9 hours  (3  hours 
after  the  maximum  land  surface  temperature 
Is  reached)  Is  shown  in  terms  of  the  two- 
dimensional  streamfunctlon  In  Fig.  3.  A 
positive  value  corresponds  to  volume  flov/ 
toward  the  land.  A rather  sharp  sea 
breeze  front  occurs  for  both  orientations 
of  the  geostrophlc  wind,  but  as  might  be 
expected  the  front  Is  sharpest  when  the 
geostrophlc  pressure  gradient  leads  to  a 
slight  opposing  wind.  In  this  case 
(ug  « -5  m/sec),  a stagnation  point  occurs 
at*^the  leading  edge  of  the  front.  The 
maximum  vertical  velocity  occurring  In 
this  front  is  about  0.09  m/sec. 
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Pig.  2.  Time  histories  at  the  shoreline 
(y  = 0)  for  Ug  ■ -5  m/sec:  (a)  V 
velocity;  (b)  U velocity. 

The  closing  of  the  circulation  cell 
occurs  over  a wide  area.  In  fact,  there 
Is  some  question  as  to  just  how  far  the 
return  flow  would  extend  out  over  the 
water  at  a later  hour  since  the  circulation 
zone  has  grown  to  such  a size  that  we  can 
no  Icnr.ci’  with  confidence  that  It  is 
unaffected  by  the  placem.ent  of  our  compu- 
tational boundaries. 
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Fig.  3.  Contour  lines  of  constant  i|»  strcamfunetlon  for  the  sea  breeze  simulation  at 
9 hrs.  after  the  land  temperature  begins  to  Increase:  (a)  u»  5 m/sec;  (b)  Ug  » -5  ra/sec. 
The  water  Is  on  the  left  (y  < 0);  the  land  on  the  right  (y  F 0). 


The  spatial  distributions  of  the  hori- 
zontal wind  occurring  at  the  same  time  as 
Fig.  3b  are  shown  In  Fig.  W.  The  sea 
breeze  extends  a little  more  than  20  km 
Inland  at  this  time.  The  return  flow  has 
a velocity  V of  about  one-half  of  the 
sea  breeze  and  is  spread  over  a layer 
about  twice  as  deep.  The  U component 
of  the  wind  near  the  leading  edge  of  the 
■ front  Is  reduced  in  amplitude  below  its 
geostrophlc  value  and  Increased  signifi- 
cantly above  its  geostrophlc  value  over 
the  water  In  the  recirculating  zone. 

Isollnes  of  turbulent  kinetic  energy 
and  humidity  mixing  ratio  are  shown  for 
this  time  and  geostrophlc  wind  direction 
In  Pig.  5.  As  expected,  the  most  turbu- 
lent region  is  in  the  sea  breeze  front, 
with  the  values  of  q*  occurring  over  the 
water  (maximum  of  about  0.l4  mVs*),  much 
lower  than  those  present  over  the  land. 

The  humidity  is  swept  along  onto  the  shore 
with  the  sea  breeze  front,  and  then  stopped 
abruptly  by  the  relatively  dry  air  over 
the  land.  If  the  computation  were  allowed 
to  run  for  several  days,  the  humidity  would 
build  up  from  the  water  surface  source  un- 
til clouds  form  at  different  points  of  the 
cycle.  Such  phenomena  must  await  future 
calculations. 

The  sea  breeze  front  at  the  time  of 
Figs.  3-5  appears  to  be  much  more  narrow 
than  can  be  well  defined  by  the  grid 
spacing  used  in  this  computation,  since 
It  extends  over  only  a few  grid  points. 

Based  on  a separate  calculation  using  a 
Ay  of  200  m over  the  reduced  domain  ex- 
tending from  200  m on  the  v/ater  side  of 
the  shoreline  to  6 km  Inland,  it  appears 
that  the  front  is  more  nearly  ]-2  km 
thick  rather  than  the  10-20  kn  implied  by 
Fig.  3b.  This  reduction  in  front  thick- 
ness also  increases  the  vertical  velocity 
within  it  by  an  order  of  magnitude.  A 
detailed,  accurate  solution  of  the  front 
structure  would  not  only  require  g.olng  to 
a finer  mesh,  but  also  relaxing  the  hydro- 
static approximation. 
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Fig.  I).  Contours  of  constant  horizontal 
velocity  at  9 hrs.  after  initialization 
for  U(T  “ -5  m/sec:  (a)  V velocity; 

(b)  U velocity. 

Figure  6 is  a comparison  of  the 
vertical  potential  temperature  Isollnes  at 
9 hours  with  those  occurring  at  18  hours. 
This  figure  demonstrates  why  there  is  such 
a large  asymmetry  between  the  day  time 
(9  hours)  sea  breeze  and  the  night  time 
(18  hours)  land  breeze  even  when  the  sur- 
face temperature  differential  is  symmetri- 
cal. During  the  day,  the  warmer  land 
temperature  leads  to  an  unstable  boundary 
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Pig*  5.  Composite  contours  of  constant 
humidity  ( — ) and  turbulent  kinetic 
energy  ( ) at  9 hrs.  after  Initiali- 

zation for  Ug  ■ -5  m/sec. 

layer  flow  which  permits  the  temperature 
differential  to  occur  over  a fairly  thick 
layer.  During  the  night,  the  surface 
cooling  Induces  a stable  boundary  layer 
flow  with  the  temperature  differential 
restricted  to  a shallow  layer.  In  previ- 
ous sea  breeze  models  occurring  In  the 
literature,  this  effect  had  to  be  simu- 
lated by  Imposing  an  empirical  eddy  dlf- 
fuslvlty  difference.  The  effect  occurs 
naturally  In  the  present  turbulent  model. 
The  cooling  over  the  water,  evident  in 
Fig.  6,  Is  a result  of  the  long-wave 
cooling  Input  with  Eq.  (2). 

^ The  horizontal  velocity  components 
(present  at  18  hours  into  the  cycle  are 
I shown  In  Pig.  7.  At  this  point,  the  land 
surface  temperature  is  280®K.  The  return 
flow  from  the  land  breeze,  l.e.,  the 
positive  V above  the  surface.  Is  almost 
nonexistent  with  a maximum  value  of 
0.35  m/sec.  A large  super-geostrophlc 
wind  parallel  to  the  coastline  exists  over 
the  breadth  of  the  land  breeze  region 
above  i/00  m.  Velocity  reversal  occurs  at  ’ 
the  surface  between  the  shoreline  and  60  km 
Inland.  It  should  be  noted  that  the  flow 
Is  not  purely  periodic  in  that  conditions 
at  t = 21/  hrs.  do  not  return  to  coincide 
with  conditions  at  t » o , 

CONCLUDING  REMARKS 

A model  of  turbulence  based  on  second- 
order  closure  has  been  applied  to  the 
atmospheric  boundary  layer  In  the  vicinity 
of  a shoreline.  The  model  is  based  on 
solving  dynamic,  partial  differential 
equations  for  the  turbulent  flux  of  mo- 
mentum, heat,  and  moisture.  These 
equations  are  obtained  by  second-order 
closure  of  the  equations  for  the  en- 
semble-averaged, single-point  moments  of 
the  fluctuating  variables.  The  numerical 
model  also  predicts  the  unsteady,  two- 
dimensional  variations  of  the  mean  values, 
and  variances  of  velocity,  potential 
temperature,  and  moisture. 


•«  Mt  •!  1*9  Ml  (— 

ItM  MMI*  Ml 

««••$  M/ltC 


Pig.  6.  Composite  contours  of  the  de- 
parture of  the  virtual  potential 
temperature  from  300“K  at  9 hrs.  (-  -) 

and  18  hrs.  ( ) after  Initialization 

of  the  simulation.  In  the  former  case, 
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Fig.  7.  Contours  of  constant  horizontal 
velocity  at  18  hrs.  after  Initialization 
for  Ujj  ■ -5  m/sec:  (a)  V velocity; 

(b)  U ^velocity. 

Calculations  have  been  made  for  a 
typical  diurnal  variation  in  the  surface 
temperature  difference  between  the  land 
and  water  with  the  geoscrophlc  pressure 
distribution  held  constant.  The  resulting 
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diurnal  variation  In  the  sea  breeze  In- 
duced by  the  strong  stability  differences 
In  the  boundary  layer  response  over  the 
land  and  that  over  the  water  produced  a 
strong  asymmetry  between  the  sea-breeze 
and  the  land-breeze  circulation  patterns. 
This  difference  was  predicted  by  the 
model  without  the  need  to  Introduce  any 
new  empirical  Information  into  the 
model.  The  sensitivity  of  the  sea  breeze 
circulation  to  the  orientation  of  the 
geostrophlc  wind  to  the  shoreline  has 
also  been  discussed. 
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